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We present a lattice Boltzmann algorithm based on an underlying free energy that allows the 
simulation of the dynamics of a multicomponent system with an arbitrary number of components. 
The thermodynamic properties, such as the chemical potential of each component and the pressure 
of the overall system, are incorporated in the model. We derived a symmetrical convection diffusion 
equation for each component as well as the Navier Stokes equation and continuity equation for 
the overall system. The algorithm was verified through simulations of binary and ternary systems. 
The equilibrium concentrations of components of binary and ternary systems simulated with our 
algorithm agree well with theoretical expectations. 



I. INTRODUCTION 

Multicomponent systems are of great theoretical and 
practical importance. An example of an important 
ternary system, that inspired the current paper, is the 
formation of polymer membranes through immersion pre- 
cipitation 1]. In this process a polymer-solvent mix- 
ture is brought in contact with a non-solvent. As the 
non-solvent diffuses into the mixture, the mixture phase- 
separates, leaving behind a complex polymer morphol- 
ogy which depends strongly on the processing conditions. 
The dependence of the morphology on the parameters of 
the system is as yet poorly understood. Preliminary lat- 
tice Boltzmann simulations of this system exist [l| . How- 
ever, this work did not recover the correct drift diffusion 
equation. A general fully consistent lattice Boltzmann al- 
gorithm with an underlying free energy to simulate mul- 
ticomponent systems is still lacking. This paper strives 
to bring us a step nearer to achieving this goal. 

There are several previous lattice Boltzmann methods 
for the simulation of multi-component systems. There 
are three main roots for these approaches. There are 
those derived from the Rothmann-Keller approach 0, Q 
that attempt to maximally phase-separate the different 
components. A second approach by Shan and Chen is 
based on mimicking the microscopic interactions [1, H, @ 
and a third approach after Swift, Orlandini and Yeo- 
mans 7, 8] is based on an underlying free energy. All of 
these have different challenges. Since we are interested in 
the thermodynamics of phase-separation we find it con- 
venient to work with a method based on a free energy. 
This allows us to easily identify the chemical potentials 
of the components. This is convenient since the gradients 
of the chemical potentials drive the phase separation as 
well as the subsequent phase-ordering. 

The challenge for the LB simulation of a multicompo- 
nent system lies in the fact that momentum conserva- 
tion is only valid for the overall system but not for each 
component separately, and diffusion occurs in the com- 
ponents. For a binary system of components A and B 
with densities p A and p B , the simulation usually traces 
the evolution of the total density p A + p B and the density 
difference p A — p B [8|] . Although this scheme is successful 
in the simulation of a binary system [{| [gj , its generaliza- 



tion for the LB simulations of systems with an arbitrary 
number of components is asymmetric. For instance, to 
simulate a ternary system of components A, B, and C 
with densities p A , p B and p c , the total density of the 
system, p A + p B + p c , should be traced, and the other 
two densities to be traced may be chosen as, e.g., p B and 
p A — p c [1J|. This approach is likely to be asymmetric 
because the three components are treated differently as is 
the case of Lamura's model [13] ■ If an LB method is not 
symmetric, it will lose generality an will only be adequate 
for special applications. In this paper, we established a 
multicomponent lattice Boltzmann method based on an 
underlying free energy that is manifestly symmetric. 



II. MACROSCOPIC EQUATIONS FOR 
MULTICOMPONENT SYSTEM 

The equation of motion for a multicomponent system 
are given by the continuity and Navier-Stokes equations 
for the overall system and a drift diffusion equation for 
each component separately. The continuity equation is 
given by 

S t p + V-J = 0, (1) 

where p is the mass density of the fluid, J is the mass flux 
which is given by J = pu, and u is the macroscopic ve- 
locity of the fluid. The Navier-Stokes equation describes 
the conservation of momentum: 

dt(pu a ) + dp(pu a up) = -dpP al 3 + dp<j al 3 + pF a , (2) 

where P a p and o a p are the pressure and viscous stress 
tensors respectively, F a is the component a of an external 
force on a unit mass in a unit volume, and the Einstein 
summation convention is used. For Newtonian fluids, the 
viscous stress tensor is given by 

<7 a( 3 = T] (dpUa + d a Uf) - ^<5 Q/ 3V • XJ^j + p B 5 a f} V • u, (3) 

where r\ is the shear viscosity, and ps is bulk viscosity, 
and d is the spacial dimension of the system. 
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Free energy, chemical potential, and pressure are key 
thermodynamic concepts to understand the phase behav- 
ior of a system. The chemical potential of each compo- 
nent can be obtained by a functional derivative as 



8T 
5n a 



(4) 



where pf is the chemical potential of component er; n a is 
the number density of component er; and T is the total 
free energy of the system. 

The pressure in a bulk phase in equilibrium is given by 



V : 



rr fir 



fl>. 



(5) 



The pressure tensor is determined by two constraints: 
Pap = P$af3 in the bulk and AP a p — ^ a n' y SJp' y every- 
where. 

In multicomponent systems, there are two mechanisms 
for mass transport: convection and diffusion. Convec- 
tion is the flow of the overall fluid, while diffusion occurs 
where the average velocities of components are different. 
The velocity of the overall fluid is a macroscopic quan- 
tity because it is conserved, but the average velocities of 
the components are not. The macroscopic velocity of the 
fluid u can be expressed in terms of the density p a and 
velocity u a of each component in the form of 



With the notation 



Au" 



(6) 



(7) 



the flux of each component can be divided into a convec- 
tion part J CTC and a diffusion part 3 ad : 



p a u a 



rad 



(8) 



Because mass conservation still holds for each compo- 
nent, the continuity equation for each component is valid: 



d t p° 



0. 



(9) 



Substituting Eq. (8J into Eq. (|9|), the convection diffusion 
equation for a component can be obtained. 



d t p a + v • r c = -v • r d . 

From Eqs. ((6|) and 0, we see that 



(10) 



(ii) 



which ensures the recovery of the continuity equation for 
the overall system. The diffusion process between two 
components is related to the difference of the chemical 
potential of the two components, which is also called the 
exchange chemical potential llj. Recognizing that the 
gradient of the exchange chemical potential determines 



the diffusion processes, we obtain a first order approxi- 
mation for the diffusion flux of one component into all 
other components as 



xad 



(12) 



where a and er' enumerate the components; pf and pf 
are the chemical potentials of components er and er'; and 
M aa is a symmetric positive definite mobility tensor. 

A simple model for the diffusion process assumes that 
a diffusion flux between two components is proportional 
to the overall density and the concentration of each com- 
ponent. Then mobility tensor can be expressed as 



M aa = k 



(13) 



where k aa is the constant diffusion coefficient between 
components a and er'. It depends on components but is 
independent of the total densities and concentration of 
each component. Substituting Eq. (fT5|) into Eq. (fT2"|) , we 
have 



V(m ct -M ct )- (14) 



Substituting Eq. p^|) into Eq. (JTUJ), the general form of 
a convection diffusion equation is obtained as 



d t p" + V(p CT u) = V V k^^-V{p- - p°'). (15) 

p 



III. LATTICE BOLTZMANN FOR 
MULTICOMPONENT SYSTEM 

To simulate a multicomponent fluid using LB we set 
up a LB equation for each component. The LBE for a 
component a of a multicomponent system is given by 



= At 



I(/-(r,i)-/f(r,t)) +J F7 



(16) 



where f?(r,t) is the particle distribution function with 
velocity Vj for component a, f ae (r^t) is its equilibrium 
distribution and F° is the forcing term of component a 
due to the mean potential field generated by the interac- 
tion of the component a with the other components. The 
main task in setting up this lattice Boltzmann method 
is to determine the correct form of the forcing term F° 
which will recover the convection diffusion equation (|15p . 

The density of each component and the total density 
are given by 



p° = £/f> 

i 

p = 5>- 



(17) 
(18) 
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The average velocity of one component a and the overall 
fluid can be defined as 



a cr 

P U Cl 



i 

pu a = ^p CT <, 



(19) 
(20) 



where u G a is the average velocity of the component cr, and 
u a is the average velocity of the overall fluid. 

The moments of equilibrium distributions for one com- 
ponent are chosen to be 

£/r = ^> 

i 

i 

S ^ j r e V la V i {i = — 8 a p + p^UaUp, 

i 

2_j f' Tev iaVipv il = — {u a 8 a p + Up8 al + v^5 a p%2\) 
i 

The moments for the forcing terms of one component are 
Y,F? = (22) 

i 

E F^ Via = p°a% (23) 

i 

J2 F °Vi a Vii3 = p»£ + ^<), (24) 

i 

^TFfViaVipVry = -p a {a'^8p 1 + ap8 ot p + a"8 a p){2h) 

i 

To utilize the analysis of the one component system we 
can establish a LB equation for the total density by defin- 
ing 

y , fi = fi ) 
a 

= F « 

a 

Y,P a < = P^ (26) 

a 

Similar to the counterparts of the one-component sys- 
tem, the moments for the overall equilibrium distribution 
function are given by 

Etf = ft 

i 

^2fiV ia = pu a , 

i 

^fiViaVi/3 = ^pS af 3 + pU a Up, 

i 

^ fiViaVipV^ = -p(u a 8p 1 + Up8 a ^ + U^8 a p) 



The moments for the overall force terms are then given 
by 

E F < = o. 

i 

i 

Using Eq. (f!M[) . we obtain 

i 

= J2(a a u p + 

cr 

= p{a a up + apu a ) + E« Au ^ + a ? A <)> (28) 



where the second term of Eq. (|28|) is of a higher order 
smallness than the first terms, and therefore does not 
enter the hydrodynamic equations to second order. For 
the third moment we have 

F l v ia v i pv il = -p{a a 8p 1 + ap8 a p + a-fSap)- (29) 



By summing Eq. (|T6|) over a, an effective LB equation 
for the total density is 



fi(r + ViAt,t + At) -fi(r,t) 



= At 



~(fi{i,t)-Mi,tj)+Fi 



(30) 



This is identical to the LB equation for a system of one 
component. Therefore, the continuity equation and the 
Navier Stokes equation of the overall fluid of a multicom- 
ponent system are recovered as 



d t p + d a (pU a ) = + O(e 3 ), 



(31) 



where U a = u a + a a At/2 is the macroscopic velocity of 
the fluid. The Navier Stokes equation for the overall fluid 
is: 

d t {pUp) + d a { P U a Up) 



-d a [ -^p8 a p 



( w 

+d a y—p(d a Up + dpU a ) 



pap 



-wd^daPUaUpUj. 



(32) 



To recover the convection diffusion equation of each com- 
ponent, we performed a Taylor expansion on the left of 
Eq. (fT6|) to second order: 



At{d t +v la d a )f° 



(Aty 



+pu a iipu- f + Q a p 7 . 



(27) 



A/ ( \{fr-f?) + F< 



(d t +v la d a ) 2 ft + 0(e 6 ) 
(33) 
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Because of the recursive nature of Eq. (|33|) . ff can be 
expressed by f° e and derivatives of f[ e as 

ft = fr+rF? 

-r(dt + v la d a )(fr + tF?) + 0(e 2 ). (34) 



Substituting Eq. J34|) into the left side of Eq. (J33J) we 
obtain 

(d t + v ia d a )(fr+TF°) 
- w(d t + v la d a ) 2 (fr + rFn + 0(e 3 ). 

= Itfr+rFf-f?). (35) 

Summing Eq. ([55)) over i gives, 

dtp" + d a {p a u a )+Td a {p a a a a ) 
'E 



u; + v la d a ) 2 (fr + tF°) = 0(e 3 ) (36) 



The first moment of ff and f? e are not identical, and the 
continuity equation cannot be obtained. Eq. (|36|) shows 
that dtp 17 + d a (p a u a ) + Td a (p a a°) is of order 0(e 2 ), and 
F? is of order O(e). Therefore d t p a +d a (p a u a ) is of order 
0(e 2 ), and we get 



>Y,(dt+vr a d a ) 2 (fr+rFn 



wd, 



dt(p a up) + d a (— + p a u a up) 



+ 0(e 3 ). 



So Eq. (f36|) can be simplified to 

d t p a + d a {p a u a ) + Td aP a al + 0{e 3 ) 



d t {p a u fj ) + d a (—+ p a u a u l3 ) 



0. (37) 



-wdfj 

Eqs. (J32]) yields 

ftCZ/j = -U a d a Up - ia Q (|<5 Q(3 ) + a/3 + 0(e 2 ). (38) 

From Eq. (35]) it follows that 

dtP° = -d a (p°U a ) + 0(e 2 ). 
Inserting Eqs. (33]) and (39]) into Eq. ([37]) we get 

p \ o 

+ p CT a/3 + 0(e 2 ). 
Substituting Eq. BD]) into Eq. (37J results in 



(39) 



(40) 



Tp a a a - Tp a a° + wd a ( — 



4* (!) 



(41) 



From this we deduce that the correct form of the forcing 
term is 



T 6 



(42) 



where the coefficient is — = 1 — =i. This coefficient ap- 
proaches 1 as At approaches 0, as one would expect from 
the continuum limit. This constitutes the main result 
of this paper. Plugging Eq. (|42]) into Eq. (|4"Tj) . we then 
obtain the convection diffusion equation 



dtP° + d a {p a U a ) = 8 a 



w 



E 



p° p ° 



d a (p a -p a ) 



43) 



The diffusion flux of component a is 



jo<l 



p"p" 



-d a (p a -P a )- 



(44) 



So that the w in Eq. (H} is equivalent to k aa ' in Eq. JT 



IV. NUMERICAL VALIDATION 

We examined the equilibrium behavior of phase sep- 
arated binary and ternary systems. We used the Flory- 
Huggins free which is a very popular model to study poly- 
mer solutions. It is given by 



+ EE^ WnCT ^ 'W> ( 45 ) 

u a' / 



where m CT is the polymerization of the component a, n a 
is its number density, and (\f is its volume fraction. It is 
defined as 



EL 
p 



(46) 



where p" is the mer density of component a and p is 
the mer density of the system, which is a constant in the 
Flory-Huggins model. To validate our algorithm we com- 
pared the binodal lines obtained by our algorithm to the 
theoretical ones obtained by minimizing the free energy. 
We used the interfacial tension parameter k = in all our 
LB simulations of binary and ternary systems, because 
there is an intrinsic surface tension in the LB simulation 
due to higher order terms , which did not appear ex- 
plicitly in the second order Taylor expansion presented 
in this paper. Since we are only evaluating the phase- 
behavior here we use a one dimensional model known as 
D1Q3. This model has the velocity set {vi} = { — 1, 0, 1}. 
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FIG. 1: (Color online) Comparison of the binodal lines for 
monomer and polymer systems show good agreement with 
theory. For m A = 1, m B — 1 we also show that the results are 
independent of the relaxation time. For the polymer system 
with m A = 10, m B — 1, the binodal points by LB with 
1/t = 0.9 the match is slightly less good for large volume 
fractions of 6 B . 
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This is an important test since all other frequently used 
higher dimensional model have this D1Q3 model as a 
projection. 

We consider two binary systems: a monomer system 
with m A = 1 and m B = 1, and a polymer system with 
m A = 10 and m B = 1. For both systems, the total den- 
sity was p — 100. Throughout this paper we choose the 
self interaction parameters to vanish: x aa — 0. The 
critical volume fractions for the monomer system are 
(f) A = 0.5 and <p B — 0.5 and for the polymer systems 
are 4> A = 0.24 and <\> B — 0.76. To induce phase separa- 
tion a small sinusoidal perturbation A<fi(x) was added 
in the initial conditions. The amplitude of the per- 
turbation is 0.1 and its wavelength is the lattice size. 
The initial volume fraction of component A is given by 
<j) A (x) = (f> A0 + A<j>(x). The initial volume fraction of 
component B is given by <fi B (x) = 4> B0 — A<p(x). 

The monomer system was simulated with different in- 
verse relaxation times. In Figure Q] we show results for 
1/t = 0.7, and 0.9. We see that the equilibrium den- 
sities have only a very slight dependence on the relax- 
ation time, although the range of stability depends no- 
ticeably on the relaxation time. The polymer system 
was simulated with only one inverse relaxation time of 
1/t = 0.9. Starting from the critical point and increas- 
ing the \ AB value for each initial condition until the sim- 
ulations were numerically unstable, we obtained a pair 
of binodal points for each initial condition. The system 
reached a stable state after about 5000 time steps. The 
measurement were taken after 50000 time steps to be sure 
that an equilibrium state had been reached. 

For the polymer system, Figure [21 shows the compar- 



FIG. 2: (Color online) The volume fraction and the chemical 
potentials of the two components for the polymer-monomer 
mixture. The parameters are m A = 10, m B = 1, \ AB ~ 0.94, 
k = 0, and 1/r = 0.9. 



ison of the total density, and the volume fractions and 
chemical potentials of each component to the correspond- 
ing theoretical values. The total density of a system in 
equilibrium by LB is essentially constant with a variation 
of Ap/p < 10~ 5 . The volume fractions of each compo- 
nent in the LB simulation agree well with the theoretical 
values. The chemical potential of each component by the 
LB simulation was very close to the theoretical value. 
The chemical potential [la, corresponding to the poly- 
mer component, varied slightly with a difference for the 
bulk values of about 2 • 10 -2 and a variation in the inter- 
face of about 4 • 10~ 2 . This is the underlying reason for 
the small deviation from the theoretically predicted con- 
centration. The potential pb was nearly constant with 
a variation of less than 10 -4 in the bulk and a variation 
of about 10~ 3 at the interface. For large values of \ AB 
this discrepancy increases leading to the noticeable vari- 
ation of the equilibrium densities of the polymer system 
as shown in Figure [TJ 

We also performed LB simulations with two ternary 



systems: a monomer system with m — 1, 



,c _ 



and ™ c 



— 1 and a polymer system with m = 10, 



and 

= 1, 



to" = 1. The x parameters for both systems were 
X AB = 3, x AC = 0-5, and X BC = 0.2. The other X pa- 
rameters were zero. The inverse relaxation time constant 
for both simulations was 1/t = 0.9. The critical point for 
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FIG. 3: (Color online) The binodal points of two ternary 
systems obtained by LB simulation. The monomer system 
has m A = 1, m B — 1, and mF = 1; the polymer system has 
m A = 10, m B = 1, and m c = 1. The parameters for both 
systems are X AB =3, x' AC = 0.5, \ BC = 0.2 and 1/r = 0.9. 
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FIG. 4: (Color online) The volume fractions and the chemical 
potentials for three component polymer system. The param- 
eters are m A = 10, m B = 1, m c = 1, x AB = 3, x AC = 0.5, 
% BC — 0.2, and 1/r = 0.9. The initial homogeneous volume 
fractions were </>a = 0.14, <j>B = 0.11, and <f>c = 0.75. 



the monomer system was (f) A ' cr = 0.32, <fr B ' cr = 0.32, and 
(j)C,cr _ 0.36. The critical point for the polymer system 
was (j) A ' cr = 0.14, (j) B ' cr = 0.11, and C ' cr = 0.75. The 
initial state of each simulation were set from the criti- 
cal points towards the end point {<j) A — 0.5, cf) B = 0.5, 
<f> = 0). Initially a small sinusoidal wave perturbation 
A(p of an amplitude of 0.1 and wavelength of the lat- 
tice size was superimposed on the initial volume fraction 
of the A component. This perturbation was subtracted 
from the B component and the C component was con- 
stant. I performed a LB simulation for each set of initial 
volume fractions and obtained the volume fractions of the 
two phases in the equilibrium state, resulting in two bin- 
odal points. The simulation reached a stable state after 
about 20,000 time steps. The measurements were taken 
after 200,000 time steps to make sure the equilibrium 
state was reached. 

Figure [3] shows the comparison of the binodal points 
by LB simulation to the theoretical binodal lines of both 
systems. The binodal points obtained by the LB simu- 
lation agree fairly well with the theoretical binodal lines 
for the monomer and polymer systems. The simulation 
becomes unstable when <p A is close to zero, i.e. when 
one component is nearly depleted. In this region the 
simulation results also deviate noticeably from the theo- 
retical binodal lines. Immediately near the critical point, 
the evolution of the system becomes extremely slow so 
the slight deviation between the binodal points obtained 
through the LB simulation and the theoretical ones prob- 
ably indicates that the LB simulation was not yet fully 
equilibrated. 

For the polymer system Figure [4] shows a comparison 
of the volume fractions and chemical potentials of each 
component. The total density of the system is again 
nearly constant with variation of less than Ap/p < 10~ 4 
in the bulk. At the interface there is a small variation 
of Ap/p < 10~ 2 . The volume fractions of each phase in 
the simulation were very close to their theoretical values. 
The chemical potential of component A was slightly dif- 
ferent in two phases with a variation of about 2 • 10~ 2 , 
while the chemical potentials of components B and C 
were much closer in the two phases with a variation of 
less than 10 ~ 4 . 



V. OUTLOOK 

We have presented a general lattice Boltzmann algo- 
rithm for systems with an arbitrary number of compo- 
nents which is based on an underlying free energy. In 
this algorithm the key thermodynamic quantities such as 
the chemical potentials of the components are immedi- 
ately accessible. It is also manifestly symmetric for all 
components. We tested the equilibrium behavior of the 
new algorithm for two and three component systems in 
each case examining both the case of monomer and poly- 
mer mixtures with an underlying Flory-Huggins free en- 
ergy. We obtained to expected phase-diagrams to good 
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accuracy and the chemical potentials were constant to 
good accuracy for the monomer systems. Polymer sys- 
tems were more challenging to simulate but still obtained 
acceptable results for m = 10. Higher polymerizations, 
however, become increasingly difficult to realize with the 
current algorithm. 

There are three directions in which we hope to extend 
this algorithm in the future. The current algorithm does 
not allow for component dependent mobility. We are 



working on developing an algorithm that can recover an 
arbitrary mobility tensor n aa . The chemical potential is 
only approximately constant. Recent progress for liquid- 
gas systems [l2| makes us hopeful that we will be able 
to ensure that the chemical potential is constant up to 
machine accuracy. And lastly we hope to extend to model 
so that it can simulate polymer systems with significantly 
larger polymerization. 
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